Isomorph invariance of Couette shear flows simulated by the SLLOD equations of 

motion 
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Non-equilibrium molecular dynamics simulations were performed to study the thermodynamic, 
structural, and dynamical properties of the single-component Lennard- Jones and the Kob- Andersen 
binary Lennard- Jones liquids. Both systems are known to be strongly correlating, i.e., have strong 
correlations between equilibrium thermal fluctuations of virial and potential energy. Such systems 
have good isomorphs, i.e., curves in the thermodynamic phase diagram along which structural, 
dynamical, and some thermodynamic quantities are invariant when expressed in reduced units. The 
SLLOD equations of motion were used to simulate Couette shear flows of the two systems. We 
show analytically that these equations are isomorph invariant provided the reduced strain rate is 
fixed along the isomorph. Since isomorph invariance is generally only approximate, a range of shear 
rates were simulated to test for the predicted invariance, covering both the linear and non-linear 
regimes. For both systems, when represented in reduced units the radial distribution function and 
the intermediate scattering function collapse for state points that are isomorphic. The strain-rate 
dependence of the viscosity, which exhibits shear thinning, is also invariant along an isomorph. Our 
results extend the isomorph concept to the non-equilibrium situation of a shear flow, in which the 
phase diagram is three dimensional because the shear rate defines the third dimension. 



I. 



INTRODUCTION 



Investigating liquids in non-equilibrium situations has 
been a matter of interest both theoretically and nu- 
merically in recent decades. Statistical mechanics pro- 
vides a link between microscopic states and equilibrium 
thermodynamics, but in non-equilibrium situations it is 
difficult to find such a linki. General formalisms for 
nonlinear response theory include the transient time- 
correlation formalism^' and the Kawasaki formalism^. 
Many theories for describing non-equilibrium liquids have 
been motivated by the theories of equilibrium situations 
and of the glassy states. Fluctuation-dissipation rela- 
tion violations^, mode-coupling theory^, and dynamical 
heterogeneity^, are examples of theories used to describe 
the behavior of systems under homogeneous shear flow. 

According to the isomorph theory^"^^ , the class of so- 
called strongly correlating liquids have isomorphs, which 
are curves in the phase diagram along which structural, 
dynamical, and some thermodynamic properties are in- 
variant when expressed in reduced units. This the- 
ory has been tested successfully experimentally^ and 
numericalljii^"— . In almost all cases studied so far 
the systems were in equilibrium, however. The above- 
mentioned interest in non-equilibrium situations moti- 
vated us to investigate whether the isomorph theory - 
or an extension thereof - holds in situations involving 
non-equilibrium steady states. To address this question 
we performed non-equilibrium molecular dynamics sim- 
ulations (NEMD) on two different systems, the single 
component Lennard- Jones (SCLJ) liquid and the Kob- 
Anderson binary Lennard- Jones (KABLJ) liquicU^"— . 
Shear flow has been previously investigated in both sys- 
tems; the SCLJ system was studied in Refs. [T9l - [2ll and 



the KABLJ system in Refs. illilil. In the NEMD sim- 
ulations reported below we focus on homogeneous flow 
generated by the SLLOD equations of motion proposed 
by Evans, Morriss and Ladd some time ago^^i^. 

The paper is organized as follows. In Sec. |ll]we briefly 
describe the theoretical basis of this work, specifically 
the SLLOD equations, and the isomorph theory. Section 
III CI proves that the SLLOD equations of motion are iso- 
morph invariant if the isomorph concept is extended to 
a three-dimensional phase diagram, where the shear rate 
is the third dimension. In Sec. IIIDI the procedure for 
generating isomorphic steady state points is explained. 
Models and simulation details are presented in Sec. IIIII 
The results of the simulations are presented in Sec. IIVI 
The paper is concluded with a summary and a brief dis- 
cussion in Sec. |Vl 



II. 



THEORETICAL BACKGROUND 



A. The SLLOD equations of motion 

Non-equilibrium molecular dynamics techniques 
(NEMD) have been used extensively to study homo- 
geneous and inhomogeneous fluids under the influence 
of different flow flelds. The case of homogeneous shear 
flow was among the first applications of NEMD^^; it 
was later generalized to elongational flows^^. Two issues 
arise when simulating shear flows. The flrst is that in a 
viscous shear flow any algorithm must satisfy the con- 
dition that the shear viscosity rj = (Jxy/'j at low strain 
rates 7 obeys the Green-Kubo linear-response relation 
(where V is the volume and T the temperature)^, 
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{Q)a.,y{t))dt . 



difference between these equations and the DOLLS equa- 
tions hes in the second term in Eqs. ([5]) and ([7]), which 
(1) has been transposed - thus the name was also "trans- 
posed" from DOLLS to SLLOD. These equations of mo- 



Here a^y is the xy element of the spatially averaged tion plus the Lees-Edwards boundary conditions-^" in con- 
junction with a Gaussian kinetic thermostat, guarantee 
that homogeneous flows in both the linear and non-linear 
regimes are generated. Excellent reviews of methods for 
simulating homogeneous flows can be found in Refs. [23l 
andlsil 

A special case of the SLLOD equations of motion is the 
Couette shear flow, where all elements of the strain-rate 
tensor are zero except one. 



stress tensor, i.e., Gxy = '^i^iFy^i/V where Xi is the x- 
coordinate of the i'th particle and Fy^i is the y-coordinate 
of the force on this particle. The second issue arising 
when simulation a shear flow is that a viscous flow nec- 
essarily generates heat. In order to simulate a steady 
viscous flow this heat must be removed, which is typi- 
cally done using a thermostat. For homogeneous NEMD 
a common method is the so-called Gaussian-- thermostat 
based on time-reversible constraint forces, which keep 
either the total energy (ergostat) or the kinetic energy 
(isokinetic thermostat) fixed. Another common method 
is the Nose- Hoover thermostat, which uses an additional 
dynamical variable to simulate the heat bath. 

Two well-known algorithms for simulating viscous 
shear flow are the DOLLS and SLLOD algorithms. The 
first homogeneous NEMD algorithm was based on the 
DOLLS Hamiltonian proposed by Hoover et al^, 



H 



DOLLS 



C/(ri,...,rAr)+^p2/2r 



E 



ri • Vv • 



(2) 

Here U is the potential energy of the system consisting of 
N particles, Vv is the gradient tensor of the macroscopic 
streaming velocity field, rui is the mass of particle i, 
and Pi are, respectively, its laboratory position and "pe- 
culiar" (thermal) momentum. The latter quantity relates 
to the velocity Ci with respect to the streaming velocity 
field according to pi = m^Ci, where Ci is given by 



V; = C; 



v(ri) 



(3) 



Via the standard Hamilton equations of motion, the 
equations generated from the DOLLS Hamiltonian are 



fi = Pi/mi -I- Ti • Vv 
Pi = Fi Vv • Pi . 



(4) 
(5) 



Here F^ is the force exerted on each particle by the sur- 
rounding particles. It was shown, however, by Evans and 
Morrissi that these equations of motion are only suit- 
able for simulating flows in the linear-response regime. 
Evans and Morriss^^ and Ladd^^ have shown that suit- 
able equations for generating flows in both the linear and 
non-linear regimes are 



i"i = Pi/rrii + ri ■ Vv 



(6) 



Pi = Fi - Pi • Vv. (7) 

The macroscopic streaming velocity fleld is assumed to 
have a linear proflle, i.e., a constant spatial gradient. The 
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(8) 



Here 7 is the shear rate, i.e., the gradient in the y- 
direction of the x-component of the streaming velocity 
fleld. Substituting the above strain rate tensor into Eqs. 
([6|) and ([7]), they take the following simple forms: 



Fj - i7Pyi . 



(9) 
(10) 



Here i is a unit vector in the positive x axis direction. 



B. The isomorph theory and its predictions 

From the instantaneous positions and momenta of all 
particles one can find the instantaneous total energy and 
pressure from 

E = if(pi,...,pAr) + [/(ri,...,r^) (11) 
pV = iVfcsT(pi,...,pAr) + T4^(ri,...,r^v). (12) 

Here if(pi, . . . , pat) and J7(ri, . . . , rjv) are the 
kinetic and potential energies, respectively, 
T(pi, . . . , Pat) is the (instantaneous) kinetic tem- 
perature given by the kinetic energy per particle, and 
W{vi,...,yn) = -l/3E»r»-Vr.C/(ri,...,r^) is the 
(instantaneous) virial, which after dividing by volume is 
the configurational contribution to the (instantaneous) 
pressure. 

Strongly correlating liquids are liquids for which a 
strong correlation exists between the constant-volume 
canonical-ensemble equilibrium fluctuations of potential 
energy U{t) and virial M^(t)2rJi,. From the fluctuations 
of potential energy and virial two parameters can be de- 
flned: 1) the density-scaling exponent 7 (the reason it is 
called the density-scaling exponent is explained in Sec. 

Ed]), 



(AWAj7) 

((AC/)2) ' 



(13) 
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and 2) the correlation coefficient R which measures how 
strongly correlating a liquid is, 



on the strain rate since this quantity defines a third state 
point coordinate. 



R 



(AWAU) 



(14) 



In both expressions the angular brackets denote an NVT 
ensemble average, AW = W{t) — {W), and AU = 
U{t) — (U). Liquids that have R > 0.9 are called strongly 
correlatingZA. For these liquids it is possible to find (ap- 
proximate) isomorphic state points in the phase diagram, 
defined as follows. Consider two state points (1) and (2) 
with densities pi and p2 and temperatures Ti and T2. 
These state points are by definition isomorphic if any 
two of their respective microscopic configurations, whose 
coordinates scale into each other according to 



^1/3^(1) _ ^1/3^(2) 
Pi '■i 



p\"vY^ (z = l,...,7V) 



(15) 



to a good approximation have proportional Boltzmann 
statistical weight factors: 



-(7(r(i',...,r5^')/A;BTi ^ (j 



126 



-U{vf\...,r^^^)/kBT^ 



(16) 



Here C12 is a proportionality constant that depends only 
on the two state points, not on the microscopic configu- 
rations. Pairs of state points in the phase diagram that 
are isomorphic fall onto the same isomorphic curve, for 
brevity termed an "isomorph" ; an isomorph is thus an 
equivalence class of state points. 

In reduced units, as a result of the proportionality 
of Boltzmann factors state points on an isomorph have 
the same dynamic, structural, and (some) thermody- 
namic quantities. Reduced units refer to the state point 
by giving lengths in units of p~^^^^ time in units of 
p^^l'^iks'T Ivn)'^!"^ , and energy in units of fc^T. The 
invariant thermodynamic quantities includ&iS the excess 
entropy (the difference between the entropy of the liq- 
uid and of the corresponding ideal gas at same density 
and temperature) and the isochoric specific heat. The 
structure is invariant along an isomorph in reduced coor- 
dinates. The equilibrium dynamic properties are also in- 
variant: Normalized time-autocorrelation functions, av- 
erage relaxation times ta = J {A{0)A{t))dt/{A'^), and 
intermediate scattering function are examples of dynam- 
ical quantities, which are invariant along an isomorph 
when expressed in reduced units^^. 

A good starting point for determining whether a quan- 
tity is isomorph invariant or not is given by Eq. (9) of 
Ref. (where the tilde signals reduced coordinates: 
r, = pi/^j..) 

{/(ri, . . . , rjv; p) = keTfiih, . • . , f^v) + 9{Q) ■ (17) 

This follows directly from the isomorph definition; Q 
labels the state point. When we generalize to non- 
equilibrium situations the term g{Q) may also depend 



C. Isomorph invariance of the SLLOD equations of 
motion 

The SLLOD equations of motion are given by Eqs. ^ 
and ([7]). In order to show that these equations are iso- 
morph invariant, one needs to substitute all quantities 
in terms of their reduced forms — the thermodynamically 
scaled dimensionless forms, denoted by tilde. Most fun- 
damental are the scaling of lengths by p~^/^ and energies 
by /cbT: thus r; = Vip~^/^ and U — UksT, whereas the 
scaling of time depends on the dynamics. If inertia is 
present, as in ordinary Newtonian dynamics and in the 
SLLOD equations, we scale the particle masses by their 
average value to; the scaling of time then follows from 
requiring invariance of the zero-strain rate equations of 
motion^''. To see this we make the above replacements in 
the SLLOD equations, denoting the time scaling factor 
as to. The scaling factor of momentum is then mp~^/^ /t^ 
(although we omit the isokinetic thermostatting term 
here for simplicity; it can also be properly expressed in 
reduced form). From the first SLLOD equation we have 
(where ifii = rrii/ {m)) 



-1/3 



Pi mp 



-1/3 



in 



■fhi torn 



-1/3 



(18) 



Clearly, all scaling factors cancel whatever choice is made 
for implying isomorph invariance independent of t^. 
From the second SLLOD equation we have 



i mp 
P 



-1/3 



-VUkeTp 



1/3 



~ _ mp 
p, Vv 



-1/3 



to' 



(19) 



where V denotes the gradient with respect to the reduced 
coordinates. Here, in order to be able to cancel the scal- 



ing factors, it is necessary that fesTp^/"^ — mp' 
or 



-1/3 



to = m/ksTp 



-1/3 



/to 



(20) 



as also shown to be the case for Newtonian dynamics in 
Ref. [l^. The reduced SLLOD equations are thus: 



= Pi/m^ -f ii ■ Vv, 



F, - p, • Vv . 



(21) 



(22) 



The isomorph invariance now follows from the fact that 
along an isomorph F^ is a unique function of reduced 
coordinates (Eq. ((IT])); this applies when the reduced 
strain rate Vv = Vv to is fixed. This means that along 
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an isomorph the strain rate must vary in such a way that 
its reduced form is fixed, see also Eq. ([50)) below. In 
the three-dimensional phase diagram parameterized by 
density, temperature, and strain rate, the isomorphs are 
one-dimensional curves. Thus it takes two parameters to 
specify an isomorph, for example the excess entropy and 
the reduced strain rate, in contrast to standard equilib- 
rium thermodynamic isomorphs, which are labelled by 
just one parameter—. 



D. Generating isomorphic state points 

It was shown recently that for strongly correlating liq- 
uids and solids, the temperature can be written as a 
product of a function of excess entropy per particle, s, 
and a function of density: T = f[s)h{p)^. Accordingly, 

one can generate curves of constant excess entropy by 

■ • s 

requiring 



-32^33 that 



Hp) 

T 



Const. 



(23) 



It is only in strongly correlating liquids that these con- 
figurational adiabats are also isomorphs, i. e., have the 
property that all the other isomorph invariants apply. 

The function h(p) is called the density-scaling function; 
its logarithmic derivative is the density-scaling exponent 
'^il^, which is also given by Eq. ([T5)l . 



7 



d\nT 
d\np 



dlnh 
din p 



(24) 



It follows that 7 only depends on the density. Note that 
in this paper 7 is always the density-scaling exponent, 
whereas 7 is always the strain rate. 

Another consequence of the isomorph theory is an 
expression for h(p) for atomic liquids with interaction 
potentials consisting of a sum of inverse power laws, 
'^('') = ^'^^ such liquids h{p) is given as 

follows^i^ 



(25) 



where the only non-zero terms are those corresponding 
to an r~" term in the pair potential. For Lennard- Jones 
liquids such as SCLJ and KABLJ, the pair potential has 
the form 



v{r) — 4e 



(7)' 



(26) 



It follows from Eq. ([25)) that the density scaling function 
can be written as h(p) = Ap'^ — Bp^; thus LJ isomorphs 
are given by an expression of the form 



Ap'^ - Bp"^ 
T 



Const . 



(27) 



Since h{p) is defined up to a multiplicative factor, 
we are free to choose a particular normalization yield- 
ing a one-parameter expression. Starting from a refer- 
ence state point (po:2o,7o), plotting the instantaneous 
virial versus the instantaneous potential energy deter- 
mines via a least-squares linear fit the value of 7 accord- 
ing to Eq. ([T3)) . denoted by 79. Following Refs. and 
I33I we can write h{p) — ap^ + (1 — ct)p^ where p = pj pQ\ 
since 7 — din h/ din p it follows that 70 = 2a + 2, i.e., 
a = 70/2 — 1. Consequently 



h = (70/2 - l)p^ - (70/2 - 2)p 



(28) 



To generate the isomorph of the reference state point we 
(repeatedly) use the equation (where T = T/Tq) 



f = h{p) , 



(29) 



and the reduced strain rate is kept fixed via (where 70 is 
the strain rate at the reference state point) 



7 = 7opi/3fi/2 



(30) 



Except for inverse-power-law systems isomorph invari- 
ance is only approximate. This means that even though 
the SLLOD equations of motion are isomorph invariant, 
it is not obvious how well isomorph invariance works in 
practice. To test this two standard strongly correlating 
systems were simulated for a range of shear rates covering 
both the linear and non-linear regimes. 



III. MODEL AND DETAILS OF SIMULATION 

The two atomic systems SCLJ and KABLJ were sim- 
ulated. In the case of the SCLJ system 500 particles in- 
teracted via the Lennard- Jones (LJ) potential Eq. ()26p . 
Reference state points for generating isomorphs are given 
by Po = 0.84, To — 0.8, with strain rates up to 2.5; all 
these quantities are given in the unit system where a = 1 
and e = 1 (LJ units). The potential was shifted and 
truncated at rcut = 3.5f7. The particles were placed in 
a cubic box, and NEMD simulations were performed us- 
ing the SLLOD equations of motion. Lees-Edwards shear 
boundary conditions were applied to eliminate effects of 
surfaces and of the small system volume. A Gaussian 
isokinetic thermostat was used to keep the temperature 
constant. The equations of motion were integrated us- 
ing the operator-splitting algorithm of Pan et alM- im- 
plemented in the GPU-accelerated MD code RUMD^; 
while RUMD, like many GPU codes, uses mainly single 
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precision floating-point arithmetic, the summation of ki- 
netic energy and similar quantities required for the isoki- 
netic thermostat was done in double precision to avoid 
unacceptable numerical drift in the kinetic energy. 

For each reference state point we plotted the instan- 
taneous virial versus the instantaneous potential energy. 
A linear regression to data according to Eq. ([T3|) gave 
7 = 5.75 and the correlation coefficient R — 0.96 at the 
zero strain rate reference state point. Within statisti- 
cal errors the value of 7 was found to be independent 
of strain rate at the reference density and temperature, 
while the correlation coefficient increases with increas- 
ing strain rate, up to about 0.99 at 7 = 2.5 (at higher 
values of 7 the system enters the so-called string phase, 
an artifact of the thermostat2§, and both R and 7 drop 
significantly). The values of R show that the system is 
strongly correlating, i.e., has good isomorphs. We stud- 
ied in detail a single isomorph obtained by increasing the 
density by 5, 10, and 15 percent with respect to po! the 
new temperatures and strain rates corresponding to each 
new density were calculated using Eqs. p9)) and (l30l) . re- 
spectively, after first determining h{p) via Eq. (1281) from 
simulations at the reference state point. 

The same procedure was applied to the KABLJ system 
with reference state points given by po = 1.2, Tq = 0.579 
(LJ units referring to the A particle parameters), and 
7 < 1.2. 1000 particles with 800 particles of type A 
and 200 particles of type B interact via Lennard-Jones 
pair potentials in a cubic box. The KABLJ potential 
parameters are as follows: uaa — 1 ctab = 0.8, ubb = 
0.88, £aa — 1, ^AB = 1-5, £bb — 0.5, and tub — riiA- We 
used the standard cut-off radius 2.5aAA- From the linear 
fit of instantaneous virial versus instantaneous potential 
energy 7 = 5.17 and R — 0.94 were determined for zero 
strain rate at the reference density. By using this value 
of 7 and Eqs. (^5)) and ([^ one can go from one state 
point to another isomorphic state point. We studied a 
single isomorph, changing density by ±5, +10, and -f 15 
percent with respect to pq. In the KABLJ system, like 
in the SCLJ system, 7 was independent of shear rate and 
R increased with increasing shear rate. 

For both systems we let the system go to its stationary 
state by running without output for some time (of order 
10^ time steps). The production phase involved of order 
lO'^-lO^ time steps, depending on the strain rate. This 
allowed for an accurate determination of structural and 
time-dependent correlation functions, as well as of the 
viscosity. 



IV. SIMULATION RESULTS 



is well known from experimental rheology of, e.g., poly- 
meric liquids^i^"— . We find the onset of the transition 
at 7 ~ 0.6 for SCLJ and 7 - 0.002 for KABLJ at the 
reference state point. 




10" 10 

Strain rate 



FIG. 1: Viscosity as a function of strain rate for (a) the SCLJ 
system at p = 0.84, T = 0.8, and (b) the KABLJ system at 
p = 1.2, T = 0.579. The transition to the nonlinear regimes 
occurs around 7 = 0.6 for SCLJ and around 7 = 0.002 for 
KABLJ. 



We used the procedure explained in Sec. IIIDI to gen- 
erate isomorphic state points for both systems, starting 
from the reference state points. Figure [5] shows the iso- 
morphic state points' density and temperature. 




FIG. 2: Density-temperature phase diagram showing four iso- 
morphic state points of the SCLJ system and five for the 
KABLJ system. The reference state points are marked with 
full symbols. 



In fluids undergoing planar Couette flow, the viscosity 
r] gives the response to the applied flow. Figure [T] shows 
viscosity versus strain rate for both SCLJ and KABLJ 
systems. At low strain rates the viscosity is constant 
(linear behavior), but at higher strain rates it starts to 
decrease. This is called shear thinning. Shear thinning 



Figure [3][a) shows the radial distribution function of 
the SCLJ system at po — 0.84, To = 0.8, for selected 
strain rates below 2. Figure [Sjb) shows the same func- 
tion for the KABLJ system at po — 1.2, Tq — 0.579, with 
strain rates below 1.2. For both systems the structure 
changes once the strain rate exceeds values correspond- 
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ing roughly to the transition to nonlinear behavior. As 
mentioned earlier, the so-called "string phases" appear at 
shear rates higher than those presented here, an artifact 
of the thermostat^. 



SCLJ 
0.84, T=0.8 



— Strain rate- 


-0 


— Strain rate- 


-0.1 


— Strain rate- 


-0.5 


— Strain rate- 


-0.9 


— Strain rate- 


-1.5 


— Strain rate: 


-2 




KABLJ 
p=1.2, T=0.579 



— Strain rate- 


-0 


— Strain rate 


-0.001 


— Strain rate- 


-0.01 


— Strain rate- 


-0.1 


— Strain rate- 


-0.5 


— Strain rate: 


=1.2 




FIG. 3: Radial distribution function g{r) of (a) the SCLJ sys- 
tem and (b) the KABLJ system at the reference state points 
with different strain rates. For clarity the radial distribution 
functions have been displaced by O.ln with n = 0, ...,5. For 
the SCLJ system there is a change of structure between strain 
rate 0.5 and 0.9, consistent with the onset of shear thinning. 
The same structure change takes place for strain rates above 
O.f for the KABLJ system, somewhat above the onset of shear 
thinning. 



Figure IDJa) shows the radial distribution functions of 
four isomorphic state points of the SCLJ system. In 
Fig-H^b) g{r) is plotted as a function of reduced distance, 
f = p^^^r. The good collapse of curves confirms the in- 
variance of structure at these strain rates that correspond 
to nonlinear flows. The same result was obtained for the 
KABLJ system; Figs. HJc) and (d) show the radial dis- 
tribution function for the five generated isomorphic state 
points in non-reduced and reduced units. 



(b) 
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— p^O.882, T=1.05, strain raie^l. 397 

— p=0.924,T=1.35, strain rate=1.607 

— p=0.965,T=1.69,strainrate=1.829 
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KABLJ 






— p=1.14,T=0.442, strain rale=0.0086 

— p=1.20,T=0.579, strain rale=0.01 

— p=1.26. T=0.741, strain rate=0.0I15 

— p=1.32, T=0.932, strain rale=0.0I31 

— p=1.38,T=1.154, strain rate=0.0I48 










I.I.I. 



2 3 
r (reduced units) 



FIG. 4: Radial distribution function for the four isomorphic 
state points of the SCLJ system in (a) non-reduced units and 
(b) reduced units, (c) and (d) Radial distribution function 
of the A particles for the five isomorphic state points of the 
KABLJ system in (c) non-reduced and (d) reduced units. To 
a good approximation the structure is invariant along the iso- 
morphs. 
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To investigate the dynamical invariance of isomorphic 
state points, we calculated the intermediate scattering 
function Fg {q, t) . In the presence of a flow, rather than 
attempting to disentangle the stochastic deviations from 
the average flow it is convenient to consider only dis- 
placements transverse to the flow direction. Following 
tradition we chose a g-value close the first peak of the 
static structure factor. Along an isomorph this q scales 
as p^^'^, and a correct comparison in reduced units must 
take this into account. We chose q = 6.81 at the reference 
density for the SCLJ system and q ~ 7.152 at the refer- 
ence density for the KABLJ system. Figure (Sja) shows 
the transverse Fs{q, t) for the SCLJ system as a function 
of ordinary time (but scaled g), while Fig.[5l^b) shows the 
same quantity as a function of reduced time. Figures[Hc) 
and (d) show the corresponding results for the KABLJ 
system. A good collapse of curves is seen when reduced 
time units are used. This shows that the dynamics are 
invariant along the isomorphs. 
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FIG. 5: Intermediate scattering function (transverse displace- 
ments) for the four isomorphic state points of the SCLJ sys- 
tem at g = 6.81(p/0.84)^'^^ as a function of (a) ordinary 
time and (b) reduced time. Intermediate scattering func- 
tion (A particles, transverse displacements) for the five iso- 
morphic steady state points of the KABLJ system at g = 
7.152(p/1.2)^'^^ as a function of (c) ordinary time and (d) 
reduced time. The collapses in (b) and (d) demonstrate iso- 
morph invariance of the dynamics in reduced units. 
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According to the isomorph theory transport quantities 
such as the reduced diffusion coefficient and the reduced 



viscosity are invariant along an isomorph. Figures EJa) 
and (c) show the viscosity along an isomorph of the SCLJ 
and KABLJ systems as functions of strain rate. By in- 
creasing the strain rate, the viscosity decreases. As men- 
tioned already this is shear thinningii22"— . In Figs. Hljb) 
and (d) the reduced viscosity rj = -q/ {p'^/^T^^'^) is plot- 
ted as a function of reduced strain rate. The isomorph 
invariance of the shear thinning phenomenon is clear. 




Strain rate 



1 

(b) 




SCLJ 






«p=0.840, T=O.S 






«p=0.882, T=1.05 






p=0.924, T=1.35 






«p=0.966, T=1.67 

















10"' 10"' 

Strain rate (reduced units) 
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10"* 10""' 10"' 10"' 

Strain rate 




FIG. 6: (a) Viscosity versus strain rate for the four isomorph 
state points of the SCLJ system; (b) reduced viscosity rj — 
v/iP^^^T^^^) versus reduced strain rate for the same state 
points, (c) Viscosity versus strain rate for the five isomorphic 
state points of the KABLJ system; (d) rj versus reduced strain 
rate for the same state points. 
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Finally, we also simulated thermodynamic quantities, 
focusing on potential energy and pressure. These quan- 
tities are not inherently isomorph invariantiiS. However, 
based on the argument* that strong correlations and the 
existence of isomorphs in non-IPL (inverse power law po- 
tential) systems is linked to a decomposition of the pair 
potential into an effective IPL part plus an almost lin- 
ear part, one expects the quantity g{Q) of Eq. (|17l) to 
depend mainly on density and negligibly on temperature 
and strain rate: The sum of all pair energies from the 
linear part of the pair potential is roughly constant at a 
given state point, and depends mainly on volume when 
different state points are considered^. The dependence 
on volume explains why quantities such as energy, free 
energy, pressure, and bulk modulus are not isomorph in- 
variant. Under the assumption that g{Q) does not de- 
pend on strain rate, the isomorph theory predicts that 
the strain-rate dependent parts of potential energy and 
virial, Uip, T, 7) - U{p, T, 0) and W{p, T, 7) - Wip, T, 0), 
are both isomorph invariant when given in reduced units. 
Notice that the same must be true for the total energy 
and pressure since the subtraction eliminates the kinetic 
terms. 




D 
D 



1 1 1 1 1 


' 1 ' ^ 


(b) 






SCLJ 




M p=0.840, T=0.8 






M.p=0.S82, T=1.05 






— p=0.924, T=1.35 






mP=0.966. T=1.67 




.•^.er-^ \ t \ t \ L 




I 



1 1.5 2 2.5 

Strain rate (reduced units) 




KABLJ 

p=1.14, T=0.442 
p=1.20,T=0.579 
— p=1.26, T=0.742 
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Figure El^a) shows the potential energy as a function 
of strain rate for the four SCLJ isomorphic state points. 
Figure El^b) plots as a function of reduced strain rate 
(7 = ^07) the strain-rate dependent part of the reduced 
potential energy, {U — Uo)/kBT, where U is the aver- 
age potential energy and Uq the average potential energy 
of the zero-strain-rate isomorphic state point. The po- 
tential energy in non-reduced and reduced units for the 
KABLJ system is plotted in Figs. (TJc) and (d), respec- 
tively. For both systems a good data collapse is seen. 



(d) 




Strain rate (reduced units) 



FIG. 7: (a) Potential energy versus strain rate for the four 
isomorphic state points of the SCLJ system; (b) The strain- 
rate dependent reduced potential energy {U — Uo) /ksT versus 
reduced strain rate to7 where (/q is the potential energy at 
zero strain rate, (c) Potential energy versus strain rate for 
the five isomorphic state points of the KABLJ system; (d) 
(U — Uo)/kBT versus reduced strain rate. 
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We did the same for pressure. Figures EJa) and (c) 
show the pressure as a function of strain rate for the 
SCLJ and KABLJ systems, respectively. Figures [HUb) 
and (d) show the corresponding strain-rate dependent 
reduced pressure (j)—po)/{pkBT) as a function of reduced 
strain rate. The cohapse is reasonable, but not as good as 
for the potential energy. We do not have an explanation 
of this. 




p=0.840, T=O.S 
p=0.882, T=1.05 - 
p=0.924, T=1.35 
p=0.966, T=1.67 _ 





FIG. 8: (a) Pressure versus strain rate for the four isomorphic 
state points of the SCLJ system; (b) the strain-rate depen- 
dent reduced pressure {jp — po) / {pknT) versus reduced strain 
rate for the same state points, (c) Pressure versus strain 
rate for five isomorphic state point of the KABLJ system; (d) 
{p — po) I (pkBT) versus reduced strain rate for the same state 
points. 
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V. DISCUSSION AND SUMMARY 

We have investigated the isoniorph theory and its pre- 
dictions for the SCLJ and KABLJ systems undergoing 
steady shear flow. Both model systems are strongly cor- 
relating (simple) liquids and known to have good iso- 
morphs in equilibrium, i.e., at zero strain rate referring to 
the standard two-dimensional thermodynamic phase di- 
agram. This paper has demonstrated that the isomorph 
concept extends to steady-state non-equilibrium situa- 
tions described by the SLLOD equations of motion, for 
which the phase diagram is three dimensional because 
the strain rate defines an extra dimension of the phase 
diagram. 

We studied structure, dynamics, and rheology in 
steady-state Couette shear flows. As expected, the struc- 
tures of both systems were unaffected by shear at very 
low strain rates, but a change of structure was observed 
as strain rate increased above the strain rate marking 
onset of nonlinear effects. The range of strain rates con- 
sidered was large enough to capture shear thinning be- 
havior. It is significant that our results include this non- 
linear regime, since the isomorph invariance of transport 
coefficients in the linear regime follows from that of the 
equilibrium properties. We obtained simulation results 
for structure studied via the pair correlation function, dy- 
namics studied via the incoherent intermediate scattering 
function, and transport quantities studied via the steady- 
state viscosity. The results show that the proposed ex- 
tension of the equilibrium isomorph theory describes well 
SLLOD steady-state non-equilibrium situations. 

In the case of the thermodynamics of the generated 
isomorphic state points, although potential energy and 
pressure are not inherently isomorph invariant, we have 
shown that the strain-dependent parts of the reduced po- 
tential energy and (to a lesser extent) pressure are invari- 
ant when considered as functions of reduced shear rate. 
It was shown by Ge et ali^ that for a dense LJ liquid 
under shear flow, the potential energy and the pressure 
can be fitted by a power-law dependence on strain rate. 



U = Uo + a-r" 



(31) 



(32) 



in which a is a common exponent that depends on density 
and temperature. They found^^ that the linear expres- 
sion a — A + BT — Cp represents well their simulations 
with A = 3.67, B = 0.69, and C = 3.35. 

The collapse seen in Fig. [7]and (to a lesser extent) and 
Fig. |8] is a consequence of the isomorph theory and the 
additional assumption that the term g{Q) does not de- 
pend on strain rate (see discussion around those figures). 



The associated master curve contains all the information 
about the strain-rate dependence of these quantities, and 
so any quantity characterizing such a master curve - for 
example a power-law exponent - is uniquely associated 
with the projection of the isomorph onto the {p, T) plane, 
which is just the equilibrium isomorph. Equivalently, 
the exponent determined by varying strain-rate at dif- 
ferent points in the (p, T) plane must be invariant along 
the equilibrium/projected isomorphs. In other words the 
theory implies that da = along an (equilibrium) iso- 
morph and, in particular, the strain-rate exponent must 
be the same for the potential energy and the pressure. 
This means that one can write 



U-Uo 



ksT 



(33) 



(34) 



in which h = a/lTf^) and b = b/{Tt^). The linear ex- 
pression of Ge et al. a — A + BT — Cp implies that a is 
constant along straight lines in the p, T plane. Accord- 
ing to the isomorph theory, however, their data would 
be even better matched by almost straight lines in the 
(Inp, InT) plane that define the isomorphs. More simu- 
lations are needed to test this prediction, but based on 
the available data we can already note the following. The 
isomorph theory implies that da = Q along an isomorph, 
i.e., BdT — Cdp = 0. Based on this one can estimate the 
density-scaling exponent from the data of Ge et al. from 
the density-temperature variation along an isomorph: 
7 = d\nT/d\np = {p/T){dT/dp) ~ (0.8/l)C/P ~ 4 
which given the uncertainties is consistent with our find- 
ings. 

We finally note that the fact that isomorph invariance 
extends beyond equilibrium situations could provide a 
powerful tool to check theories of non-equilibrium be- 
havior. This is because isomorph invariance imposes a 
constraint on the temperature and density dependence 
of transport coefficients, and any general theory for these 
must result in an isomorph invariant expression for the 
reduced transport coefficients. This would be analogous 
to the "isomorph filter" for theories of the dynamics of 
viscous liquids approaching the glass transitiorkiS.. 
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